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Abstract 

Previous simulation and experimental studies have shown that some grain boundaries 
(GBs) can couple to applied shear stresses and be moved by them, producing shear defor- 
mation of the lattice traversed by their motion. While this coupling effect has been well 
confirmed for symmetrical tilt GBs, little is known about the coupling ability of asymmet- 
rical boundaries. In this work we apply a combination of molecular dynamics and phase 
field crystal simulations to investigate stress-driven motion of asymmetrical GBs between 
cubic crystals over the entire range of inclination angles. Our main findings are that the cou- 
pling effect exists for most of the asymmetrical GBs and that the coupling factor exhibits a 
non-trivial dependence on both the misorientation and inclination angles. This dependence 
is characterized by a discontinuous change of sign of the coupling factor, which reflects a 
transition between two different coupling modes over a narrow range of angles. Importantly, 
the magnitude of the coupling factor becomes large or divergent within this transition re- 
gion, thereby giving rise to a sliding-like behavior. Our results are interpreted in terms of a 
diagram presenting the domains of existence of the two coupling modes and the transition 
region between them in the plane of misorientation and inclination angles. The simulations 
reveal some of the dislocation mechanisms responsible for the motion of asymmetrical tilt 
GBs. The results of this study compare favorably with existing experimental measurements 
and provide a theoretical ground for the design of future experiments. 
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1 Introduction 



Recent research has led to the recognition that many grain boundaries (GBs) in crystalline materi- 
als can couple to applied shear stresses and are moved by them in a manner similar to dislocation 
glide Q]-|9]|. During this "coupled" GB motion, the boundary produces shear deformation of 
the lattice and causes relative translation of the grains parallel to the GB plane. The coupling 
effect is characterized by a coupling factor (3 defined as the ratio of the tangential grain transla- 
tion velocity vu to the normal GB velocity v n : (3 = Vu/v n . Coupling is considered perfect if (3 
is a geometrically determined number depending only on crystallographic characteristics of the 
boundary and not on its velocity or the driving stress. 

Initially observed experimentally for low-angle GBs in Zn ifTOl [TT| . the coupling effect has 
recently been demonstrated for many high-angle GBs in a number of metallic and non-metallic 
materials [|6]-|9l [T2l - l20l . The coupling can be responsible for the stress-induced grain growth 
in nano-crystalline materials Iff6l I2TT - I241 and can influence the nucleation of new grains during 
recrystallization 11251 . Molecular dynamics (MD) simulations have been effective in providing 
insights in the atomic mechanisms, geometric rules and dynamics of coupled GB motion [|2]- 
l4l l2Tll26] - l3Ti . It has also been found that the phase field crystal (PFC) methodology developed 
over the recent years I132I4451I is well capable of reproducing equilibrium and non-equilibrium GB 
properties. In particular, it predicts reasonable values of GB energies 11421 and has been able to 
reproduce the phenomenon of GB premelting 1138113911 as well as non-trivial structural transitions 
at GBs at high homologous temperatures [44J. In addition, it has been shown to reproduce the 
coupling effect for both two-dimensional (2D) [|43ll and 3D symmetric tilt boundaries 051. 

Most of the experiments as well as simulations conducted so far have been focused on sym- 
metrical tilt GBs. The less-studied case of asymmetrical tilt boundaries is more complex but 
poses new interesting questions. For example, low-angle symmetrical tilt GBs are known to 
move by collective dislocation glide in parallel slip planes 11461 . By contrast, low-angle asym- 
metrical GBs are composed of at least two different types of dislocations. One would expect that 
the dislocations gliding in intersecting slip planes could block each other and prevent the coupled 
motion. In fact, the impossibility of coupled motion of asymmetrical boundaries was suggested 
in the classical paper by Read and Shockley [|46l . Yet, recent MD simulations B3U 071 1481 . 
bicrystal experiments (6l |9l and the observation of coupled GB motion in poly crystalline mate- 
rials suggest that this may not be the case. Even less is known about geometric rules of coupling 
or migration mechanisms of large-angle asymmetrical GBs. 

Based on purely geometric considerations, it was suggested that deviations of the GB plane 
from symmetric inclinations should preserve the coupled motion for both low and high-angle 
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GBs [4J. Furthermore, the coupling factor was predicted to be independent of the inclination 
angle as long as the coupling mode remains the same (see discussion of the multiplicity of cou- 
pling modes below) @|. The MD simulations [47] of [001] tilt GBs in copper with the tilt angle 
9 = 18.92° and several different inclinations showed that (3 did vary with the inclination angle 
within 10 to 20%. However, the boundary remained coupled over the entire angular range of 
inclinations as predicted. On the other hand, for large-angle boundaries with 9 = 36.87°, two 
symmetrical and two asymmetrical boundaries were tested and the coupling factors were found 
to be different in both magnitude and sign (two /3's were positive and two negative) [47]. Zhang 
et al. [31J studied stress-driven motion of asymmetrical [001] tilt GBs in Ni with 9 = 36.87° and 
9 = 33.36°. They observed coupled motion (with varying sign of (3) in some cases and sliding 
in other cases. The inconclusive and often conflicting results of the previous studies point to the 
need for a more detailed and systematic analysis of coupling of asymmetrical GBs. 

In this paper we report on simulations of stress-driven motion of a large set of asymmetrical 
[001] tilt GBs by applying two complementary methodologies: MD and PFC. The MD simula- 
tions are conducted for specific materials (FCC copper and aluminum) and provide quantitative 
information about the mechanical stresses needed for driving the GB motion at different temper- 
atures. The MD simulations are also well suited for studying atomic-level mechanisms of GB 
migration by examining atomic trajectories. A weakness of the MD approach is that the time 
scale is limited to tens of nanoseconds, preventing access to diffusion-controlled processes such 
as dislocation climb. 

The PFC methodology overcomes, in principle, the latter limitation by permitting simulations 
on diffusive time scales, thereby describing both dislocation glide and climb [|34l . This gives us 
the opportunity to get a glimpse into the possible GB behavior in the long-time regime which 
is currently unaccessible by MD simulations. However, one limitation of the PFC method is 
that the number of peaks of crystal density waves does not directly correspond to the number of 
atoms, which is not generally conserved. While this permits a description of dislocation climb 
phenomenologically, it remains unclear how to define and control the vacancy concentration 
in the PFC model. Theoretical attempts have been made to incorporate vacancies explicitly 
in the PFC model, but they do not fully resolve this issue since climb remains possible even 
for a vanishing vacancy concentration [|49l . As a consequence, even though the geometrical 
aspects of the glide-mediated conservative GB motion are well modeled by PFC quantitatively, 
the description of non-conservative climb-mediated motion remains largely qualitative. 

Despite this limitation, we find here that the MD and PFC approaches reveal remarkably sim- 
ilar coupling behaviors as a function of the misorientation and inclination angles, including the 
existence of a narrow transition region between two coupling modes and a sliding-like behavior 
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in this transition region. Hence, a comparison of the MD and PFC simulations helps us shed 
light on aspects of GB behavior that do not depend sensitively on detailed atomic mechanisms. 
Finally, in order to be able to explore efficiently the entire parameter space of the misorientation 
and inclination angles, we restrict our PFC simulations to 2D GBs between square lattices. We 
make use of the PFC model of Ref. BHJ, which energetically favors the square and FCC lattices 
in 2D and 3D, respectively. Even though the 2D square geometry is simpler than the 3D FCC 
structure studied in the MD simulations, we are able to quantitatively relate the Burgers vector 
character of the GBs in the PFC and MD simulations owing to the fact that the square lattice 
rotated by 45° is identical to one plane of the FCC lattice. Hence, the choice of the 2D geometry 
is not highly restrictive for the purposes of the present study. 

2 Simulation geometry 

The geometry of asymmetrical [001] tilt GBs between FCC crystals is illustrated in Fig. |TJ Be- 
cause the orientations of the GB plane and the tilt axis are fixed, the system has two geometric 
degrees of freedom. These degrees of freedom can be associated with the tilt angle 9 and the 
inclination angle 0. The tilt angle is defined as the misorientation angle between the [001] di- 
rections in the grains, with 9 = corresponding to a single crystal and 9 > to the [001] i axis 
rotated in the counter-clockwise direction relative to [001]2. The inclination angle is defined 
as the angle between the GB plane and the internal bisector between the [001] directions in the 
grains. We take > if the bisector is rotated in the counter-clockwise direction relative to the 
GB plane. The case of = with 9^0 corresponds to a symmetrical tilt GB. 

Due to the fourfold symmetry of the FCC lattice, all distinct GB structures can be found in the 
angular domain {0 < < 7r/4, < < 7r/4}. However, we are interested in not only the GB 
structures but also their orientations relative to the laboratory coordinate system xyz (Fig.[T]). In 
order to include all possible orientations of the GB structures, we consider an expanded domain 
{0 < 9 < tt/2, -tt/4 < < vr/4}. 

Symmetry analysis gives the following orientation relationships between the GB structures: 

(0,0)^(0,-0), (1) 



flL-r 

^ (9O°-0,45°-0), </>> 0, 



(90° - 9, -45° - 0) , < 0. 



(2) 
(3) 



4 



Eq. ([T]) states that the GB structures with the same 9 but opposite signs of <p are mirror reflections 
of each other across the GB plane (x, z). Eqs. ^ and ^ state that replacement of the angles 9 
and <p by the complementary angles (90° — 9) and (±45° — 0) reflects the GB structure across 
the mirror plane m x normal to the x-axis. In particular, the inclination angles = ±45° for a 
given 9 produce identical symmetrical tilt boundaries which are m x -reflections of the symmetri- 
cal boundary with = and the tilt angle (90° — 6). As will be discussed later, these symmetry 
relations impose certain conditions on the response of the GBs to applied shear stresses. The 
symmetry relations ([l])-([3]) also apply to the 2D GBs modeled by the PFC method with an appro- 
priate choice of the lattice unit cell. 

3 Multiplicity of coupling modes 

When a GB executes coupled motion, its dislocation content shears the lattice swept by its motion 
and simultaneously rotates the lattice to align it with the lattice of the growing grain. This 
combination of shear and rotation produces relative translation of the two grains parallel to the 
GB plane by the amount f3L, where L is the normal GB displacement. As was discussed in 
previous work [|2l-Bl, the coupling factor (3 is a multi- valued function of the crystallographic 
angles of the boundary. There are two ways to understand the origin the multiplicity of coupling 
factor. 

One way is to recognize that the shear deformation produced by the GB depends on its dislo- 
cation content, which, for a general GB, is defined by the Frank-Bilby equation [50]. The latter is 
known to have multiple solutions due to the crystal symmetry, leading to the multiplicity of dis- 
location content of the GB. Different dislocation contents produce different deformations of the 
lattice, which is manifested in different coupling factors observed during the boundary motion. 
This results in the existence of multiple coupling modes of the same GB, each corresponding to 
a different solution of the Frank-Bilby equation. 

Another interpretation of the coupled motion focuses on the lattice rotation step. Consider 
a tilt GB as an example. To ensure continuity of the lattice of the growing grain, the receding 
lattice must be rotated around the tilt axis by the angle ±9 (the sign depends on the direction of 
GB motion). However, if the lattice possesses n-fold rotation symmetry around the tilt axis, then 
rotations by the angles ±9 + (2ir/n), ±9 + 2(2n/n), etc., also produce physically identical states 
of the growing grain. But these different rotation angles lead to different relative translations of 
the grains and thus different coupling factors. 

For the [001] tilt axis in a cubic material, the fourfold symmetry generates four possible 
coupling modes, with the coupling factors (3 = 2 tan (9/2 + irk/ A), k = 0, 1, 2, 3. In reality, 
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only two of them, corresponding to the smallest magnitude of (3, have been observed in MD 
simulations [|2]-|4l and experiments Il6l [T5l [T81 - I201 . These two modes are referred to as (100) and 
(110) type and are characterized by the coupling factors 



respectively. Note that these two coupling factors have different signs and describe GB motion in 
opposite directions in response to the same shear. For symmetrical tilt boundaries, the coupling 
factors obtained by simulations EJ-H1 and experiments ll6l lT5l[T&l - l20l accurately follow the (100) 
mode for angles 9 < 36° and the (110) mode for angles 9 > 36°. At the critical angle of approx- 
imately 36°, (3 abruptly changes from one mode to the other and can exhibit a "dual behavior" 
in which the boundary switches back and forth between the two modes [EJ-HJ. It should be noted 
that this switching angle is not prescribed by any symmetry requirements and its exact value may 
depend on the material and/or the crystal structure. In particular, in a previous MD study of Cu 
flU, an examination of the gamma- surfaces for different coupling modes (Fig. 24 in [4]) revealed 
that the slip responsible for the (100) mode is more difficult than the slip corresponding to the 
(110) mode. The lower Peierls-Nabbaro (PN) barrier for slip in the (110) mode is consistent with 
the fact that this mode spans a larger range of misorientation, and therefore that the switching 
angle is less than 45°. 

We now discuss the coupling factors of asymmetrical tilt GBs. Suppose, as suggested in 
[ErHl, the coupling factor for a given mode is independent of the inclination angle. Then one 
can hypothesize that the expected angle dependence of (3 would look as shown schematically 
in Fig. [2] Only the (100) and (110) modes are considered and their coupling factors are repre- 
sented by two surfaces. The cut between the surfaces corresponds to the discontinuous transition 
between the modes accompanied by a reversal of sign. The exact shape of this cut is unknown 
and it is drawn in this figure with the only requirement that it respect the symmetry relations 
Q-Q. When applying these relations, it was taken into account that the reflection of the GB 
structure across its plane (x, z) does not affect the coupling factor. By contrast, reflection of the 
GB structure across the mirror plane m x normal to the x-axis reverses the sign of (3. In particular, 
if for symmetrical boundaries with = the mode switching occurs at a tilt angle of 9 = 36°, 
then for symmetrical boundaries with <p = ±45° it occurs at 9 = 90° — 36° = 54°. Accordingly, 
all boundaries with 9 < 36° are expected to have the same positive coupling factor, given by 
Eq. Q, regardless of 0. Likewise, all boundaries with 9 > 54° are expected to have the same 




(4) 



and 




(5) 
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negative coupling factor given by Eq. ^ regardless of (p. For boundaries with 36° < 9 < 54°, 
(3 is positive at and near = but is expected to switch to negative values when approaching 
= 45° and <p = —45°. Although the exact switching angles may depend on the material, this 
analysis gives rather definitive qualitative predictions that can be tested by simulations. 

It should be emphasized again that the shape of the cut between the surfaces in the putative 
diagram of Fig.[2j and even its existence, are at this point hypothetical. It is the goal of this paper 
to test this diagram by simulations. Our simulation results reported later in this paper will reveal 
that this diagram is qualitatively correct in identifying the overall shape of regions of opposite 
signs of /3. However, we will see that it fails to capture the fact that the magnitude of f3 depends 
on the inclination angle and can become very large near the discontinuous boundary between the 
two coupling modes. 

4 Methodology of atomistic simulations 

Atomic interactions were modeled using embedded-atom method (EAM) potentials fit to exper- 
imental and first-principle data for Cu and Al |[5Tl l52ll . Both potentials accurately reproduce 
physical properties of these metals that are important in the context of this study. In particular, 
they predict accurate values of the elastic constants and stacking fault energies as well as dislo- 
cation core structures. The melting points predicted by these potentials are T m = 1327 K for Cu 
and T m = 1040 K for Al (the experimental values are 1357 K and 933 K, respectively). It should 
be noted that Al and Cu have significantly different stacking-fault energies and thus different 
dislocation core splittings and mobilities. The inclusion of both metals under the same method- 
ology was intended to lend this work more generality and facilitate comparison with possible 
future experiments. 

The GBs were created by constructing two separate crystals and joining them along a plane 
normal to the y-direction (Fig. [T]). Periodic boundary conditions were imposed in the x- and z- 
directions parallel to the GB plane. To achieve commensurability of lattice plane periodicities in 
the grains in the x- and ^-directions, the angles 9 and <fr were chosen so that to create a coincident 
site lattice (CSL) and align one of the CSL planes with the GB plane. The dimensions of the 
simulation block in the x- and ^-directions comprised integral numbers of CSL periods, which 
completely eliminated coherency strains. In the y-direction, the grains were terminated at free 
surfaces. Several atomic layers near each surface were exempt from the MD process and were 
used only to control the boundary conditions. Namely, atoms in surface layer 2 were fixed in 
their perfect lattice positions during all simulations, whereas atoms in surface layer 1 were fixed 
only in the y- and ^-directions and were displaced with the same constant velocity in the x- 
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direction. All other atoms were dynamic. Approximate dimensions of the simulation block were 
L x ~ 100-200A, L Y « 370A and L z « 37A and the total number of atoms was 1 x 10 5 to 
2 x 10 5 . 

The equilibrium GB structure was obtained by the energy minimization procedure described 
in ll53l . Prior to the MD simulation, the block was uniformly expanded by the thermal expansion 
factor at the chosen simulation temperature. This expansion was intended to eliminate thermal 
stresses inside the grains. The thermal expansion factors for EAM Cu and Al were known from 
separate calculations [|27l I5T1 . The MD simulations were performed in the canonical (NVT) 
ensemble with temperature controlled by a Nose-Hoover thermostat. A 2 femtosecond time 
integration step was used throughout this study. After temperature reached the target value, the 
GB was equilibrated by an isothermal anneal for a few hundred picoseconds. The equilibration 
was followed by a production run in which the surface layer 1 was moved parallel to the re- 
direction with the speed of v = 0.5 m/s, imposing a shear strain rate of about 10 7 s _1 . Shears 
in both positive and negative x-directions were implemented for each boundary. The production 
runs took up to 50 nanoseconds but were often terminated earlier when the GB reached one of 
the surface layers (Fig.[TJ). Multiple snapshots of the simulation block containing the coordinates 
of all atoms, their energies and other relevant information were saved during the simulations. 

The following technique was applied for tracking the GB motion. First, an orientation param- 
eter, A, was assigned to each individual atom in a given snapshot. To this end, vector positions 
Tij of n ~ 50 nearest neighbors of a chosen atom i were compared with positions f m (a) of 12 
first- nearest neighbors of an atom in the perfect FCC lattice rotated by an angle a around [001]: 

n 12 3 
3=1 m=l k=l 

Here a is the equilibrium lattice constant and index k runs over three Cartesian components 
of a vector. The quantity ipi(a) was calculated for two angles a corresponding to the chosen 
rotations of the grains relative to the coordinate system. The angle with the larger ipi{a) provided 
a better match between the actual and ideal lattice orientations and was assigned to atom i as its 
orientation parameter A. Fig. [3] illustrates this procedure by showing the atoms assigned to the 
grains by the bright and dark colors. 

After partitioning the atoms between the grains, the mean GB position h was calculated as 
h = L Y N/N , where N is the number of atoms in grain 1 and iVo is the total number of atoms 
in the simulation block. The plots of h versus time were used for calculation of the coupling 
factor (3. The first and last 15 A of the plot were disregarded to eliminate the effects of the elastic 
strain and interaction with the surface layers. The remaining part of the plot was linearized by a 



{Tijk ~ fmk{0t)Y 



(6) 
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least-mean-square fit and the slope S of the regression line was used to compute (3 = v/S. 

Various visualization methods were applied for examining the structures of moving GBs. 
In particular, for GBs composed of discrete dislocations their slip traces contained information 
about the dislocation movements and reactions. The deformation fields method proposed in ll54ll 
was implemented and applied for visualization of dislocations and their slip traces during the GB 
motion. 



5 Methodology of PFC simulations 

The PFC model is generally formulated as an evolution equation for the dimensionless crystal 
density field 4>(r,t) defined as the departure of the atomic number density n(r,t) from some 
reference value n normalized by that value, ^(r, t) = (n(r, t) — no) /no- Here PFC simulations 
are carried out in the grand canonical ensemble with a constant chemical potential /i, so that ij) is 
a non-conserved field. We use the dynamical formulation [|40ll 

where the second order partial derivative with respect to time, d 2 ?p/dt 2 , has the advantage of 
rapidly relaxing the elastic field via propagation of wave-like modes that mimic phonons in a 
real solid. The first order partial derivative with respect to time damps those modes and also 
models diffusive processes. The total free-energy of the system is expressed in terms of the 
functional 

F = [ drf + F ext , (8) 



where / is the free-energy density of the bulk crystal chosen to have the form [HTTl 

/ = \ H + (V 2 + l) 2 [(V 2 + Qlf] ) + £ (9) 

and F ext specified below models an "external" potential used to shear the system. The form of 
/ couples two sets of crystal density waves with different reciprocal lattice vectors, where Q 1 is 
the ratio of the magnitudes of those vectors. This form models FCC lattices in 3D by coupling 
[111] and [200] reciprocal lattice vectors with Qi = a/4/3, and square lattices in 2D by coupling 
[10] and [11] reciprocal lattice vectors for the choice Q 1 = \J2, which is adopted here. 

As in the MD simulations, the GBs were created by constructing two separate crystals and 
joining them along a plane normal to the y-direction (Fig. [T]). We use identical definitions and 
sign conventions for the tilt and inclination angles as in the MD simulations. Periodic boundary 
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conditions were applied in the x-direction, with the angles 6 and and the dimension Lx chosen 
to be equal to an integral number of CSL periods in order to eliminate coherency strains in the 
crystals. The dimension Ly was chosen to be large enough to obtain GB behaviors that are 
independent of Ly, with typical values in the range of 100 lattice spacings. 

To enable comparison between the PFC and MD results, we replaced the primitive square unit 
cell of the 2D PFC lattice with edge a by an expanded unit cell rotated by 45° with edge y/2a, 
which contains 2 atoms (Fig. [4]). As a result, the PFC structure becomes identical up to a scaling 
factor to one of the (002) crystal planes of the 3D FCC structure. With this choice of the unit cell, 
the Burgers vectors of the two types of dislocations present in the MD and PFC lattices have the 
same Miller indices: (100) and 1/2(110) in MD and (10) and 1/2(11) in the PFC simulations, 
respectively. Henceforth, all crystallographic indices related to the PFC structures will be given 
with respect to the expanded unit cell. Fig. [4] illustrates 1/2(11) dislocations forming a low- 
angle symmetrical tilt GB. The Burgers vector is given by the closure failure of the Burgers 
circuit drawn around one of the dislocations. 

Since the PFC method does not model a solid-vacuum interface, we cannot use free surfaces 
as in the MD simulations to shear the bicrystal. However, we can mimic free surfaces by choosing 
the chemical potential to vary spatially over narrow strips near the top and bottom surfaces of 
the crystal for the geometry depicted in Fig. [T] This is accomplished by choosing the chemical 
potential /i to vary spatially along the direction y normal to the GB from a value ji s which favors 
the solid phase in most of the bulk of the sample to a value [L\ which favors the liquid phase near 
the top and bottom surfaces. With y varying from to Ly, where corresponds to the bottom 
surface, the choice 

A* = W + ^ ~ ^ (tanh [(y — Ly + &)/£] - tanh [(y - &)/£]) (10) 

melts the top and bottom surfaces a distance b into the sample if the parameter £ is chosen 
much smaller than b. To determine the values of Hi and /i s , we first compute the solid- liquid 
phase diagram by a standard common tangent construction where the free-energy of the solid 
phase is computed by expanding the crystal density field in terms of the two sets of [10] and 
[11] density waves as described in PTTl . This construction yields the value of the equilibrium 
chemical potential, [i E , for solid-liquid coexistence. The solid (liquid) phase is favored for jj, > 
fiE < He) on the side of the phase diagram where [ie is negative and the average solid density 
is larger than the liquid density. We choose fi s = 0.9/ig and /ij = 1.1/ze for all simulations 
reported here. 



10 



A shear is imposed by choosing 

Fext = J dxdy [G(y - d) (ip(x,y,t) - ip (x ~ vt,y)f 

+ G(y -L y + d) (i/i(x, y, t) - ^(x + vt, y)f] , (11) 

where G(y) = exp(— y 2 /2a 2 ) / \phxo 2 is a normalized Gaussian of width a and ipo(x, y) corre- 
sponds to the analytical expression for the equilibrium ip field for a perfect single crystal approxi- 
mated as a superposition of the [10] and [11] sets of crystal density waves (Eq. (87) in IHTTO . Since 
the dynamics tends to drive the system towards a minimum of the grand potential F — [i J drijj, 
the terms proportional to (ip(x, y, t) — ipo{x ± vt, y)) 2 in the integrand of F ext tend to entrain ip 
in a way that is equivalent to pulling the crystal in opposite directions along x at velocity ±v. The 
Gaussian kernels enforce that this pulling only takes place in narrow strips of width a centered 
at a distance d into the sample from its bottom and top surfaces. In general, d needs to be chosen 
much smaller than L y and slightly larger than the width b of the surface melted layers so that the 
pulling is applied to solid strips near the bottom and top surfaces that are not melted. 

All the simulations were carried out with a = 0.1 and v = 1.25 x 10~ 4 , where 2v is the 
relative velocity of the two crystals (cf. Eq. <[TT[)). A series of simulations of symmetrical tilt 
boundaries and all the simulations of asymmetrical tilt boundaries were carried out with e = 0.25 
for which \i E = —1.33215 (recall that pL s = 0.9/x.e and pLi = 1.1/i^inEq. ([TO])). We also repeated 
the series of simulations for symmetrical tilt boundaries with e = 0.05 and thus he = —0.593103. 
This choice was motivated by the fact that the discontinuous transition between the two different 
coupling modes is expected from the previous MD work flU reviewed above to occur at a critical 
misorientation that is controlled by the relative magnitudes of the PN barriers to dislocation 
motion in different coupling modes. For the higher e value (e = 0.25), the PN barrier turns out to 
be substantially larger for the (10) than 1/2(11) dislocations 115511 . Since low-angle symmetrical 
tilt boundaries with small 9 are composed of (10) dislocations, the coupled motion is expected 
to switch from the (10) to the (11) coupling mode at a relatively low misorientation, i.e., at a 
misorientation that is large enough for (10) dislocation cores to overlap but smaller than 45°. By 
contrast, for the lower e value, the PN barriers are much smaller and of comparable magnitude 
for the (10) and 1/2(11) dislocations. In this case, one would expect the transition between the 
two coupling modes to occur at a misorientation angle close to 45°. Our PFC simulation results 
reported below will confirm these expectations by showing that the switching angle between 
the coupling modes is close to 20° when the PN barrier is much larger for (10) than for 1/2(11) 
dislocations (e = 0.25) and comes close to 45° when both barriers are negligibly small (e = 0.05). 

We note that the calculation of the PN barriers in PFC simulations requires choosing the 
shearing velocity v sufficiently small for the viscous-like stress associated with grain translation 
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(absent in the MD simulations and real solids) to be much smaller than the Peierls stress. Compu- 
tation of this stress and reformulation of the PFC model to eliminate this spurious viscous stress 
will be discussed elsewhere ll55Tl . What is important for the present study is that v was chosen 
small enough for the results not to depend on this PFC artifact. 

The dynamical equations were solved using a pseudo-spectral scheme with a mesh spacing 
Ax = Ay = 27r/8 and a time step At = 0.1. This scheme extends the one described in 
11391 to treat the second order partial derivative in time in Eq. ([7]). Details of this scheme will 
be described elsewhere P3i Since the edge of the square unit cell with one atom per cell is 
a = 2tc in dimensionless PFC units (and \pla for the expanded unit cell with 2 atoms used for 
comparisons with MD), this choice of mesh spacing corresponds to having 8 x 8 = 64 mesh 
points per unit cell of the square lattice. For the parameters used to melt and pull the top and 
bottom surfaces of the crystal, we chose b = 5a = lOn, £ = 2a, a = a = 2n, and d = 12a. 

6 Equilibrium grain boundary structures 

Examples of equilibrium GB structures studied in this work are given in Figures [5] and [6] Fig.^a) 
shows an asymmetrical boundary with relatively small angles of 9 = 16.26° and = 14.04°. 
In the structure obtained by the MD simulations, individual GB dislocations can be easily dis- 
tinguished and their Burgers vectors can be readily determined by the standard Burgers circuit 
construction. This boundary has a periodic structure comprising six (100) dislocations and four 
1/2(110) dislocations in each period. The core of each dislocation can be considered as a stack 
of identical structural units representing capped trigonal prisms in three dimensions. Their 2D 
projections appear as kites and are outlined in the figure in green color. The PFC simulations 
give a similar structure of this boundary, with the same number of dislocations of each type 
(Fig. [5jb)). The exact positions of the dislocations are different and vary with temperature and 
time. However, the dislocation content predicted by both methods is identical. In fact, all GBs 
with 9 = 16.26° studied in this work were found to be mixtures of these two types of dislocations 
in proportions dictated by the inclination angle. These proportions were found to be precisely the 
same in both MD and PFC simulations. In the particular case of = 0, the boundary becomes 
symmetrical and its structure represents an array of (100) dislocations. Likewise, both methods 
confirm that the symmetrical boundaries with <p = ±45° are composed of 1/2(110) dislocations. 

As examples of high-angle GBs, Figs.|6ja,b) show the structures of symmetrical boundaries 
with 9 = 36.87°. This lattice misorientation produces a CSL known as E5, £ being the recip- 
rocal density of coincident sites ll50ll . The structures shown in this figure are for the inclination 
angles = and = 45°, corresponding to the GB planes (310) and (210), respectively. Both 
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structures are formed by topologically identical kite- shape structural units similar to those in 
Fig. [5ja). Such units are separated by one atomic bond when = and are connected head-to- 
tail when = 45°. The rows of such structural units running parallel to the tilt axis are similar 
to dislocation cores in low-angle boundaries and can be also considered as GB dislocations in 
high-angle boundaries. 

Figure[6]^c) illustrates a typical structure of an asymmetrical £5 GB, with the inclination angle 
of = 14.04°. The kite-shape structural units "clustered" together or separated by a bond can 
be considered as patches (facets) of the symmetrical £5 (210) (0 = 45°) and £5 (310) (0 = 0) 
boundaries, respectively. 

The PFC simulations show similar structural trends of high-angle GBs. However, a detailed 
one-to-one comparison between the PFC and MD structures is limited because of their different 
dimensionality. In particular, the characteristic kite-shape structural units forming the 3D GBs 
include sites located in adjacent (002) planes. Such structural units have no analog in the 2D 
PFC structures containing only one layer. 

7 Mechanisms of grain boundary motion 
7.1 MD results 

The mechanisms of GB motion will be discussed for relatively low-angle boundaries in which 
the evolution of individual dislocations could be reliably traced. Examination of MD snapshots 
revealed that the GB motion at low temperatures (0.3T m to 0.5T m ) was accomplished by dislo- 
cation glide assisted by dislocation rearrangements and reactions. As mentioned in Section [TJ 
stress-driven motion of symmetrical tilt GBs occurs by collective glide of identical dislocations 
along parallel slip planes (3l |4|. This motion does not require dislocation rearrangements or re- 
actions. This mechanism was indeed observed in our simulations as illustrated in Fig. |7fa) for 
the boundary with 9 = 16.26° and = 0. 

In the case of asymmetrical GBs, the two types of dislocations forming the GB structure 
belong to different slip systems. After a period of time, the dislocations gliding in intersecting 
slip planes can create locked configurations preventing further GB motion. Nevertheless, our 
MD simulations have shown that the dislocations usually find a way to glide past each other 
without completely locking themselves. Fig.^b) illustrates the simultaneous glide of (100) and 
1/2(110) dislocation arrays during the motion of an asymmetrical boundary with 9 = 16.26° 
and = 38.66°. Note how the dislocation slip traces cross each other multiple times without 
locking. 
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Two mechanisms were identified by which the dislocations could avoid blocking each other 
while preserving the total Burgers vector. These mechanisms involve dislocation reactions and 
dislocation avoidance, respectively. 

The dislocation reaction mechanism is similar to the one observed in the recent MD study of 
shrinkage and rotation of isolated cylindrical grains 11281 . Fig. [^schematically illustrates a typical 
dislocation reaction process in which a single 1/2 (110) dislocation propagates through an array 
of (100) dislocations. The structure shown in this figure represents a typical asymmetrical GB 
with relatively small angles 9 and <p < 9. At each step of this process, the (100) dislocation on the 
immediate right of the 1/2 (110) dislocation dissociates in two 1/2 (110) 's. One of the product 
dislocations glides over a short distance (comparable to the dislocation spacing) and recombines 
with the initial 1/2 (110) dislocation to form a new (100). The remaining product dislocation is 
similar to the initial 1/2 (110) but is located a short distance to its right (compare configurations 
(a) and (g) in Fig. [8]). This new configuration continues to propagate to the right by repeating the 
same steps: reaction of the 1/2 (110) dislocation with a neighboring (100) and recombination 
with one of its products. This process looks as if the 1/2 (110) dislocation migrated along the 
boundary to the right (compare the initial (a) and final (h) configurations in Fig. [8]). However, it 
is only the 1/2 (110) Burgers vector that propagates over many steps, whereas each individual 
dislocation glides only over a short distance. Note that this dislocation propagation process 
produces a slight displacement of the dislocation array down. Multiple dislocation passes can 
produce significant GB displacements. The remarkable feature of this mechanism is that it does 
not require dislocation climb, despite the fact that the propagating dislocation has a Burgers 
vector component normal to the GB plane. 

A similar chain or reactions can propagate a (100) dislocation through an array of 1/2 (110) 's 
or (100)'s, or a 1/2 (110) dislocation through an array of 1/2 (110)'s. Such dislocation reactions 
provide a mechanism for rapid redistribution of dislocation content over the GB without altering 
the total dislocation content or relying on diffusion-controlled mechanisms such as climb. These 
dislocation reactions prevent the formation of locks and simultaneously provide a mechanism for 
GB motion. 

The second mechanism was the dislocation avoidance. When the ratio of the numbers of 
the two type of dislocations was large, we observed that the minority dislocations tended to 
lag behind the majority dislocations and then return to the boundary when a suitable gap was 
available. This process is illustrated schematically in Fig. [9] for the dislocation ratio 2:1. First, 
the majority dislocations move forward while the minority dislocations are left behind. This 
separation of dislocations creates gaps in the array of majority dislocations. (Such gaps represent 
elastically distorted perfect lattice regions between terminations of dislocation arrays and can be 
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considered as disclination dipoles [56J). Then, when the slip planes of the minority dislocations 
come to alignment with the gaps, they quickly glide forward and fill the gaps, recreating the initial 
GB structure. As evident from the geometry of this process (Fig. [9]), each minority dislocation 
fills a gap left by a neighboring minority dislocation. Thus, the net result of this process is 
similar to the dislocation reaction mechanism: both processes lead to relative translations of two 



dislocation arrays parallel to the GB plane without interfering with each other. Fig. 10 illustrates 
the dislocation avoidance mechanism for a particular GB with 9 = 16.26° and 4> = 2.73°. This 
Figure shows that the minority dislocations do not necessarily fill all gaps simultaneously but can 
fill them one or several at a time. 

The chains of dislocation reactions represent the dominant mechanism responsible for the 
motion of asymmetrical GBs. The dislocation avoidance mechanism was observed only for 
GBs with small inclinations and inclinations close to ±45°, which contained large ratios of the 
different dislocation types (e.g., 10:1). In a small number of simulations at temperatures below 
0.3T m , the dislocations formed strongly locked configurations. Such locks eventually triggered 
generation of new dislocations which initiated plastic deformation of the grains. 

Similar mechanisms were found to operate in high-angle GBs, with rows of kite-shape struc- 
tural units playing the role of dislocation cores. 

7.2 PFC results 

In the PFC simulations, the dislocations forming the GBs with relatively low angles could be 
readily identified and followed during the GB motion. As expected, symmetrical boundaries 
with = and <p = ±45° migrated by glide of identical dislocations in their respective slip 



planes. Fig. 11 shows a typical dislocation trace during the motion of a GB with 9 = 16.26° and 



45°. This boundary is composed of 1/2 (11) dislocations which glide in (11) planes. A more 



detailed picture of this process is illustrated in Fig. 12 The dislocation moves by a conservative 
process in which the structural units forming its core are continually distorted and converted to 
perfect-lattice units left behind the dislocation. 

For asymmetrical GBs, however, the migration mechanisms are more complex. Typically, 
the majority dislocations glide in their respective slip planes as before, whereas the minority 
dislocations progress parallel to the same slip planes as the majority. This process is illustrated 
in Fig.[T3]by dislocation traces in the asymmetrical GB with 9 = 16.26° and = 30.2°. This 
particular boundary contains one (10) dislocation per every three 1/2 (11) dislocations. The 
majority dislocations move by perfect glide along (11) planes containing their Burgers vector. 
The minority (10) dislocations move parallel to the same (11) planes, a process which cannot 
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be explained by perfect glide. Indeed, the dislocation Burgers vector has a component normal 
to (11) planes, suggesting that the dislocation motion must be accompanied by some amount of 
climb. 

The latter conclusion was verified by observations of the detailed core structure of the minor- 
ity dislocations during their motion. As illustrated in Fig.[l4j the number of sites (atomic density 
peaks) in the dislocation core is not conserved. In this particular example the core loses one site, 
whereas in other cases it could gain sites. As noted in Section|5} PFC simulations model an open 
system, in which the sites are not conserved locally or globally. The PFC simulations represent 
the material's behavior at high temperatures approaching the melting point. Therefore, the con- 
tinual creation, disappearance and redistribution of the sites can be interpreted as occurring by 
diffusion of vacancies. This interpretation can explain the motion of the minority dislocations 
along the majority slip planes by a combination of glide and climb. It should be emphasized that 
this mechanism was not, and could not be observed by the MD simulations because of the short 
time scale. 

As mentioned earlier, the 2D character of the PFC methodology precludes a direction com- 
parison of the structural units in high-angle GBs with their MD counterparts. Furthermore, the 
non-conservative nature of the PFC simulations often obscures unambiguous interpretation of 
atomic movements in complex GB structures. These complications prevented us from a more de- 
tailed PFC study of structural evolution and migration mechanisms of asymmetrical high-angle 
GBs. 



8 Orientation and temperature dependencies of the driving 
stress 

Previous MD simulations ll26l 1271 have shown that the resistance of symmetrical tilt GBs to 
coupled motion is relatively small and is due primarily to the stick-slip friction associated with 
nucleation of disconnection loops. The present MD simulations indicate that for asymmetric 
GBs, the resistance to motion is much greater and is caused by the need to avoid or overcome 
locked configurations between the dislocations gliding in intersecting slip planes. This resis- 
tance can be characterized by the average shear stress a xy required for sustaining a constant GB 
velocity. 



Typical time dependencies of the shear stress are shown in Fig. 15 for an asymmetrical GB 
in Al. The initial rise of the stress is due to accumulation of elastic deformation until the stress 
reaches a level sufficient for sustaining the GB motion. The values of a xy reported below were 
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obtained by averaging over the steady state portion of the simulation run. At high temperatures 
the stress behavior is uniformly noisy, whereas low temperatures reveal multiple peaks charac- 
teristic of lock-unlock dynamics. It should be mentioned that in the case of perfectly regular 
stick-slip behavior, the average stress depends on the system size Ly in the direction normal the 
GB plane ll27ll . However, the size-dependent correction to the stress decreases as 1/Ly and for 
the large system sizes studied here is small. 



The angle dependence of the average shear stress is plotted in Fig. 16 for a series of GBs 
with the same tilt angle 9 = 16.26° but varying inclination angle. Observe that the stress is 
highest for the "most asymmetrical" GBs with the inclination angles |0| ~ 20°-30°. By contrast, 
the stresses required for moving the perfectly symmetrical GBs arising at <p = and ±45° are 
almost an order of magnitude smaller. 



Fig. 16 demonstrates that the driving stress strongly decreases with temperature while pre- 
serving the same trend of the angle dependence. A more detailed temperature dependence of the 
stress is illustrated in Fig. [TTJfor an asymmetrical GB in Al. The rapid decrease of the stress 
with temperature reflects the thermally activated nature of the mechanisms responsible for the 
dislocation reactions, avoidance and unblocking. At high temperatures approaching the melting 
point, the driving stress is a factor of ten smaller than at room temperature. 



9 Orientation and temperature dependencies of the coupling 
factor 

For symmetrical GBs, the coupling factor depends on the tilt angle 9 and temperature. As dis- 
cussed in Section [3j previous MD simulations and experiments have revealed that (3 is a multi- 
valued function of 9 and exhibits a discontinuous transition between two coupling modes. This 
transition was also confirmed in the present MD simulations (not shown here). The existence of 
two coupling modes and a discontinuity were also reproduced by the PFC simulations as illus- 



trated in Fig. 18 The jump of the coupling factor occurs at an angle close to 20° for e = 0.25, 
which is the value of e used in all the simulations of asymmetrical GBs, and close to 45° for 
e = 0.05. This difference reflects the e-dependence of the PN barriers of the two types of dis- 
locations as discussed earlier. Note the excellent agreement with predictions of the geometrical 
model of coupling 01, in which the two branches of the plot are described by Eqs. Q and ([5]). 
For asymmetrical GBs, we first discuss the 9 = 16.26° GBs which were studied in greatest 



detail. Fig. 19 reports the MD results for such boundaries in Cu and Al, showing two temper- 



atures in each case. As indicated in Section [3j for 9 below the critical angle of approximately 



17 



36°, (3 is expected to remain a positive constant equal to the ideal value given by Eq. ([4]). In- 
stead, the simulations show that (3 of asymmetrical GBs varies with the inclination angle and, to 
a lesser degree, with temperature. The deviations from the ideal value of (3 are positive for some 
inclinations and negative for others. In Al the deviations can be as high as a factor of two. Al 
is dominated by positive deviations whereas in Cu both positive and negative deviations occur to 
nearly equal extent. 

Similar results were obtained by PFC simulations (Fig. [20]). In this case, the angle of discon- 



tinuity of symmetrical GBs is about 20° (Fig. 18), thus for 9 = 16.26° the coupling factor was 
again expected to remain constant. Instead, it varies with in a manner reminiscent of that in Al 
(cf. Fig. [T9fb)). The deviations from the ideal (3 are always positive and reach a peak at about 
±30°. It should be emphasized that, despite the significant deviations of (3 from its ideal value 
in both MD and PFC simulations, the coupling factor remains positive. The positive sign of (3 
indicates that the coupling mode remains the same at all inclination angles, which is consistent 
with our geometric analysis in Section [3j 

Fig. |2T|a) shows the MD results for high-angle GBs with 9 = 36.87°. In this case, the 
coupling factor is expected to be negative at = and change sign to positive values as \(j>\ 
increases towards 45° (Section [3]). This behavior is indeed confirmed by the simulations, despite 
the significant scatter of the points and deviations from the ideal values of (3 (Fig. |2"T|a)). The 
scatter strongly increases in the discontinuity regions, in which the magnitude of (3 becomes very 
large. Importantly, the PFC simulations reveal a very similar behavior of the coupling factor 
(Fig.[2ljb)). The PFC plot is smoother than the MD plot and clearly reveals a "divergence" of (3 
{(3 — > ±oo) in the narrow region where the sign changes. 

Due to the high computational efficiency of the PFC method, it was possible to perform 
a large set of simulations with various misorientation and inclination angles. The results are 
summarized in Fig. [22] as a diagram in the coordinates 9-<\> showing positive and negative values 
of (3 by different symbols. The boundary between the positive and negative regions has been 
drawn in this figure by hand and is intended to be a trend line. Despite the obvious scatter of the 
points, the shape of this line is in qualitative agreement with the geometric predictions shown 
in Fig. [2] as far as the overall shape of the two regions of positive and negative (3 is concerned. 
In particular, this diagram clearly shows that the discontinuity between the two coupling modes 
known from the previous work |2]44l [6l [331 [T8l - l20l exists not only for symmetrical boundaries 
but also extends over the full range of inclination angles. However, sections of this diagram at 
a constant misorientation angle in the MD (Figs. 19 and [2T|a)) and PFC (Figs. 20 and [2T|b)) 
simulations clearly show that the magnitude of (3 is strongly dependent on inclination in a way 
that cannot be predicted by the geometrical considerations which were used to construct Fig.[2j 
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10 Discussion 



The goal of this paper was to investigate the effect of the inclination angle on stress-driven 
motion of asymmetrical tilt GBs. To ensure generality of the results, the MD simulations were 
conducted over the full range of inclination angles, a wide temperature range, and in two different 
metals. Furthermore, the same GBs were also studied by PFC simulations, a recently developed 
methodology which provides access to atomic-level processes on diffusive time scales Il32]4l31l . 

The MD and PFC methods are complementary to each other and their combination offers 
an efficient approach to multiscale problems such as GB motion. MD simulations are suitable 
for studying the effects of temperature, strain rate and stress on GB dynamics Il26l 1271 while 
simultaneously providing detailed information about atomic mechanisms. However, dislocation 
climb and other diffusion-controlled processes are beyond the time scale of MD simulations 
(tens of nanoseconds). The PFC approach is less quantitative and cannot be used for tracking 
all details of atomic movements. However, the material is modeled in the long-time regime in 
which atomic diffusion can readily occur. This gives access to GB migration mechanisms (such 
as dislocation climb) that would otherwise not be seen. 

This study has shown that at a fixed misorientation angle 9, the coupling factor (3 varies 
with the inclination angle <\> (Figs. 19} 20} [21) . This observation does not confirm the geometric 
prediction 01 that (3 is a function of 9 only. But there are other theoretical predictions that were 
tested and verified by this work. The most important of them is that most GBs in materials are 
coupled and, unless subject to imposed constraints, can be moved by applied stresses [ffl. Indeed, 
the overwhelming majority of asymmetrical GBs tested here coupled to shear stresses and were 
moved by them. This finding clearly demonstrates that the coupling effect is not an attribute of 
only symmetrical GBs which were predominantly studied in previous work. 

Furthermore, our simulations show that the two coupling modes found previously in [001] 
symmetrical tilt GBs [|2Tl4l [61 [T51 [T81 - I201 continue to exist in asymmetrical GBs. The discontin- 
uous transition between the modes also continues to exist for asymmetrical GBs as illustrated in 
Figs. [2]and 22 In the vicinity of the transition between the modes, the magnitude of the coupling 
factor becomes very large if not divergent (Fig. 21 ). In other words, a GB caught between the two 
coupling modes responds to applied shears by a processes which appears like sliding. Further 
investigations are needed to determine whether this response is a manifestation of a true sliding 
process or of the "dual behavior" with alternation between the two coupling modes as suggested 
on p. 4974 of [ffl. 

The long-standing mystery [46] related to asymmetrical boundaries is how the dislocations 
gliding in intersecting slip planes avoid blocking each other (Fig. [7]). Strong interactions between 
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such dislocations can explain the high stresses needed for the coupled motion of asymmetrical 



GBs (Fig. [16]). At the same temperature and GB velocity, the shear stress required for steady- 
state motion of an asymmetrical GB can be an order of magnitude greater than for symmetrical 
boundaries. Several mechanisms have been found by which the dislocations alleviate the locks, 
most notably the dislocation reactions (Fig. [8]) and dislocation avoidance (Figs. [9] and [10]). When 
diffusion is allowed, the locks can also be overcome by dislocation climb as suggested by the 



PFC simulations (Section 7.2). The operation of these complex mechanisms responsible for the 
unlocking of the dislocations can explain the deviations of the coupling factor from the ideal 



geometrical value (Figs. [19} [20] and [2T]). Such deviations are predominantly positive, suggesting 
the existence of a sliding component along with coupling. However, a better understanding of 
the origin of the deviations requires further studies. 

One might think that the deviations of the coupling factors from the ideal values could be 
caused by an additional driving force arising due to the elastic anisotropy of the lattice. For 
asymmetrical GBs, the elastic anisotropy creates a difference in elastic strain energies in the 
grains which may drive the GB motion. In fact, this driving force underlines one of the methods 
for studying GB migration by atomistic computer simulations, see e.g. [5] for review. However, 
because this driving force is quadratic in stress, large stresses need to be created in the grains in 
order to drive GB motion over a non-negligible distance. The stresses used in our study were 
smaller and unlikely to explain the above deviations. This is also evident from Fig. [17] showing 
the temperature dependencies of the coupling factor and stress. At temperatures above 500 K, f3 
does not change within the scatter of the data points and remains higher than the ideal value. In 
the same temperature interval, the stress drops by about an order of magnitude, which rules out 
the possibility that the deviation from the ideal (3 was caused by additional GB motion induced 
by the stress. 

It is interesting to compare the predictions of this modeling study with experimental data. Our 
main finding that most of asymmetrical GBs are moved by applied shear stresses is well consis- 
tent with experimental observations of stress-driven grain growth in nano-crystalline materials 
H6l|23l|24l|57l. Coupled motion of asymmetrical tilt GBs was also observed in experiments on 
bicrystalline samples subject to applied shear loads ttHHl. In particular, Molodov et al. |6l studied 
stress-driven motion of [001] tilt GBs in Al bicrystals with different crystallographic parameters. 
One of the experiments was performed on a bicrystal with 9 = 17.4° and <\> = 19.1°. The bound- 
ary was found to couple to the applied stress (see Figure 7 in O) and move with a coupling 
factor of 0.39, which is higher than the ideal geometrical value 0.31 computed from Eq. Q. To 
our knowledge, this is the only published value of coupling factor for an asymmetrical tilt GB in 
Al. This value can be directly compared with our simulations. Indeed, from the MD simulation 
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results for Al shown in Fig. [T9~fo), the boundary closest to the experimental conditions is one 
with 9 = 16.26° and 4> = 18.44°. The computed coupling factor for this boundary is 0.40 at the 
temperature of 500 K and 0.38 at 900 K. Although the experimental temperature corresponding 
to /3 = 0.39 was not specified in [6], our prediction is in excellent agreement with experiment at 
both temperatures. Our simulations indicate that the overestimated coupling factor is caused by 
the dislocation reactions as discussed earlier in this paper. It should also be pointed out that the 
overestimated (relative to the geometric prediction) (3 value found in (6]| is well consistent with 
the general trend found in this work by both MD (Fig.[T9)^b)) and PFC (Fig. 20) simulations. 



The symmetry analysis (Sections|2]and[3]) and the simulations reported in this paper provide a 
theoretical basis for designing future experiments. In particular, for [001] tilt GBs in FCC metals, 
the following experiments would provide the most direct test of our theoretical predictions. First, 
for a given material, a series of symmetrical tilt GBs should be studied over the entire range of 
misorientation to determine the critical switching angle 9 C between the different coupling modes 
(100) and (110). This series of experiments would parallel the experimental study of symmetric 
tilt GBs in Al where 9 C was estimated to be between 30.5° and 36.5° [fT9l . After this critical 
angle has been established, the following series of experiments should be able to reveal the three 
distinct regimes predicted in this paper: 

1. Measurements of (3 at a fixed tilt angle 9 < 9 C and varying inclination angle 0. The 
coupling factor is expected to remain positive but may deviate from the value predicted by 
Eq. ([4]) (most probably in the positive direction). 

2. Measurements of (3 at a fixed tilt angle 9 > (90° — 9 C ) and varying inclination angle 4>. 
The coupling factor is expected to remain native but may deviate from the value predicted 
by Eq. ([5]) (most probably in the negative direction). 

3. Measurements of (3 at a fixed tilt angle 9 C < 9 < (90° — 9 C ) and varying inclination angle 
0. The coupling factor is expected to switch from positive at small <p to negative as 4> 
approaches ±45°. At angles close to the sign change, the boundary may exhibit a sliding - 
like behavior. 

It is important to emphasize that the critical switching angle between the two coupling modes is 
essential in determining the range of misorientation, 9 C < 9 < (90° — 9 C ), which is predicted to 
exhibit the most interesting and novel dependence of the coupling on inclination. As indicated 
earlier, the specific critical angle of 9 C = 36° mentioned above refer to EAM Cu and can be 
different for other metals. 
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11 Conclusions 



In conclusions, we have shown by computer simulations that the coupling effect exists for the vast 
majority of symmetrical as well as asymmetrical tilt GBs. Furthermore, as for symmetrical GBs, 
the coupling of asymmetrical boundaries can occur in multiple models with coupling factors 
(3 that can have different signs. The magnitude of f3 is generally strongly dependent on the 
inclination angle in a way that cannot be fully predicted from purely geometrical considerations. 
The most dramatic manifestation of this dependence is the sharp increase in the magnitude of 
(3 in a range of angles near the boundary between regions of opposite signs of (3. This confers 
GBs with a sliding-like behavior that could potentially have a strong influence on mechanical 
properties of poly crystalline materials. 

Furthermore, we have found that the motion of asymmetrical GBs can be mediated by several 
different processes. In the MD simulations where vacancy diffusion was negligible, the motion 
of GB dislocations on different slip planes was accommodated by dislocation reactions and/or 
avoidance. In contrast, in the PFC simulations which incorporate diffusive processes generally 
occurring at high temperatures, dislocation climb facilitated collective motion of dislocations 
with different Burgers vectors, allowing the GB to avoid locks and move smoothly. 

Importantly, the dependence of the coupling factor (3 on GB bicrystallography was found 
to be strikingly similar in the MD and PFC simulations, despite the mentioned differences in 
the atomistic details of GB migration. Both simulation methods predict similar shapes of the 
regions of opposite signs of (3 in the parameter space of angles. These shapes are consistent 
with considerations of crystal symmetry and the different PN barriers of the dislocations moving 
in different slip planes. Both methods predict diverging magnitudes of (3 near the boundary 
separating the regions of positive and negative values of (3. We therefore expect these basic 
features of asymmetrical GBs to pertain to a wide range of materials, temperatures and other 
physical conditions. 

Finally, our simulation results are in encouraging agreement with experiments on stress- 
driven GB motion in poly crystalline materials and bicrystalline samples. Furthermore, they are in 
excellent agreement with the recently measured experimental value of the coupling factor for an 
asymmetrical tilt GB in an Al bicrystal JH. When discussing the problem of coupled motion of 
asymmetrical boundaries, Molodov et al. [6] remark: "Further investigations, especially molec- 
ular dynamics simulations are obviously needed to clarify the mechanisms of this phenomenon, 
specific atomic rearrangements, dislocation processes, and reactions involved in the process of 
boundary migration." We hope that the present study has made a step in this direction and will 
motivate further experiments. 
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Figure 1 : Geometry of the MD simulation block employed in this study. 
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Figure 2: Schematic dependence of the coupling factor f3 on the tilt angle 6 and the inclination 
angle </>. (a) Surface plot with two sheets representing the (100) and (110) modes of coupling. 
The coupling factor is discontinuous along the cut. (b) Projection of the surfaces on the 6*-0 
plane, showing the areas of positive and negative values of (3. 



29 




Figure 3: Illustration of partitioning of atoms between the two grains using the orientation pa- 
rameter method. The image was taken from MD simulations of the Cu GB with 9 = 16.26° and 
4> = 14.04° at 500 K. The bright and dark colors designate atoms assigned to different grains. 
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Figure 4: PFC simulated structure of the symmetrical tilt GB with 6 = 80°. The GB plane 
is horizontal and its approximate position is indicated by an arrow. The structure is composed 
of 1/2(11) dislocations whose cores are marked by yellow triangles. A Burgers circuit drawn 
around one of the dislocations is shown in white color. The unit cell of the 2D lattice and a 
perfect (11) crystal plane passing between the dislocations are indicated. 
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Figure 5: Equilibrium structure of the asymmetrical GB with 6 = 16.26° and <p = 14.04°. (a) 
Obtained by MD simulations of Cu and Al. The open and filled circles mark atomic positions 
in alternate (002) planes. The structural units forming the dislocation cores are outlined, (b) 
Obtained by PFC simulations. The dislocation cores are marked by yellow triangles. Both struc- 
tures are composed of six (100) dislocations and four 1/2(110) dislocations in each structural 
period. 
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Figure 6: Examples of equilibrium structures of £5 GBs (9 = 36.87°) obtained by MD sim- 
ulations of Cu and Al. (a) Symmetrical (310) boundary with = 0. (b) Symmetrical (210) 
boundary with <\> = 45°. (c) Asymmetrical boundary with <fi = 14.04°. The kite-shape structural 
units of the GB structures are outlined. Selected crystal planes are labeled for clarity. 
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Figure 7: Traces of gliding dislocations in MD simulations of coupled GB motion in Cu at 500 
K. The dislocation traces are revealed using the micro-rotation vector method from Il54l. The 
green and blue colors represent different localized lattice deformations. The initial and final GB 
positions are indicated, (a) Symmetrical tilt GB with 9 = 16.26° and = moves by collective 
glide of (100) dislocations, (b) Asymmetrical tilt GB with 9 = 16.26° and (p = 38.06° moves by 
collective glide and reactions of (100) (minority) and 1/2(110) (majority) dislocations. 
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(a) Initial state 



(b) Dissociation and glide 



(c) Recombination 



(d) Intermediate State 



(e) Dissociation and glide 



(f) Recombination 



(g) Intermediate State (h) Final state (i) Summary of propagation 

Figure 8: Schematic of dislocation propagation along a GB by a chain of dislocation reactions. 
The dislocation notation is the same as in Figure[5} The dislocation pairs undergoing dissociation 
and recombination reactions are encircled, (a) Initial state, (d) and (g) Intermediate states, (h) 
Final state, (i) Summary of the whole process. 
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Figure 9: Schematic illustration of the dislocation avoidance mechanism for an asymmetrical 
GB moving upward. The dislocation notation is the same as in Figure [5j (a) Initial dislocation 
structure comprising (100) (majority) and 1/2(110) (minority) dislocations. The dashed lines 
indicate the dislocation slip planes ahead of their motion, (b) The majority dislocations move 
forward while the minority are left behind, creating gaps in the GB structure, (c) When the gaps 
are aligned with the slip planes of the minority dislocations, the latter move forward and fill the 
gaps, recreating the initial GB structure in a new position. 
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Figure 10: MD snapshots illustrating the dislocation avoidance mechanism for an asymmetrical 
GB in Cu with 9 = 16.26° and <\> = 2.73°. The GB moves upward, (a) Initial GB structure, (b) 
The minority dislocations 1/2(110) are left behind creating gaps between the majority disloca- 
tions (100). (c) When proper alignment is reached, two of the minority dislocations catch up with 
the boundary and fill the gaps. The remaining minority dislocations subsequently fill other gaps 
(not shown). The dislocations are visualized by constructing small Burgers loops around their 
cores and color-coding a selected projection of the closure failure. The simulation temperature 
is 500 K. 
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Figure 11: A typical trace of a 1/2(11) dislocation in PFC simulations of coupled motion of the 
symmetrical GB with 9 = 16.26° and <p = 45°. The dashed lines indicate (11) and (10) crystal 
planes in the advancing grain. The dislocation moves parallel to (11) planes consistent with the 
glide mechanism. The distances are in Angstroms. 
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Figure 12: PFC observation of 1/2(11) dislocation glide during coupled motion of the symmet- 
rical GB with 9 = 16.26° and = 45°. Frames (a), (b) and (c) show sequential configurations 
of the boundary moving down, with the current GB position shown by an arrow. The dislocation 
core is marked by a yellow triangle. The white line connects a selected set of sites (the same in 
all three frames) located in a (11) plane. One of the sites is labeled by an open circle for tracking. 
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Figure 13: Typical dislocation traces in PFC simulations of coupled motion of an asymmetrical 
GB with 9 = 16.26° and <\> = 30.3°. The dashed lines indicate (11) and (10) crystal planes in the 
advancing grain. Both (10) and 1/2(11) dislocations move parallel to (11) planes, suggesting 
that the motion of (10) dislocations involves climb. The distances are in Angstroms. 
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Figure 14: PFC observation of (10) dislocation motion during coupled migration of an asym- 
metrical GB with 9 = 16.26° and <\> = 30.3°. Frames (a) and (b) show sequential configurations 
of the boundary moving down, with the current GB position shown by an arrow. The dislocation 
core region is outlined by five selected sites marked by open circles and connected by yellow 
lines. In (a), the yellow polygon encircles a sixth site, whereas in (b) this site disappears, provid- 
ing evidence of dislocation climb. The white lines outline some of the (20) planes as a guide to 
the eye. 
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Figure 15: Typical time dependencies of the shear stress during coupled motion of an asymmet- 
rical GB in Al with 6 = 16.26° and = —18.44° at three different temperatures. 
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Figure 16: Steady- state shear stress as a function of inclination angle for GBs with the tilt angle 
9 = 16.26°. (a) Results for Cu. (b) Results for Al. The simulation temperatures are indicated in 
the legends. 
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Figure 17: Shear stress (a) and coupling factor (b) as functions of temperature for an Al GB with 
9 = 16.26° and = —18.44°. In (b), the horizontal line indicates the ideal coupling factor. 
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Figure 18: Coupling factor of symmetrical tilt GBs as a function of misorientation angle obtained 
by PFC simulations with e = 0.25 and e = 0.05. The lines indicate theoretical predictions H| 
for two coupling modes. 
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Figure 19: MD simulation results for the coupling factor of asymmetrical GBs in (a) Cu and (b) 
Al with 6 = 16.26° as a function of the inclination angle (p. The simulation temperatures are 
indicated in the legends. The horizontal line indicates the ideal coupling factor. 
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Figure 20: PFC results for the coupling factor of asymmetrical GBs with 9 = 16.26° as a function 
of the inclination angle 0. The horizontal line indicates the ideal coupling factor. 
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Figure 21: MD simulation results for Cu at 800 K (a) and PFC simulation results (b) for the 
coupling factor of asymmetrical GBs with 9 = 36.87° as a function of the inclination angle 0. 
The horizontal lines indicate the ideal coupling factors for two coupling modes. The shaded 
stripes indicate approximate regions in which the coupling factor switches between the modes. 
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Figure 22: Summary of PFC calculations of the coupling factor for asymmetrical GBs. The 
diamond and circle symbols indicate positive and negative (3, respectively. The line outlines the 
approximate boundary between the two coupling modes with different signs of (3. 
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